In [1]:

    

import numpy 
import matplotlib
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt


%matplotlib inline
from matplotlib import rcParams
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16


    

In [2]:

    

nx=30 # the number of nodes in lattice direction 
dt=1
dx=1
x1=0.
x2=1.0 # the length 
alpha=1.22e-3 # heat diffusion coefficient 
D=1 # the dimension of the problem 
mstep=70000 # the number of time step
omega=1./(0.5+(alpha*D))     #Chapman-Enskog  dt=1 and dx=1 
Tleft=100.0    # left wall temperature

k=3 # k=0,1,2    <==== c(2)===c(0)====> c(1)


    

In [3]:

    

x=numpy.linspace(x1,x2,nx+1)
w=numpy.zeros(k)      # witghting factor
T=numpy.zeros(nx+1)   # Temperature matrix
f= numpy.zeros((k, nx+1))   # distribution function
feq= numpy.zeros((k, nx+1)) # Equilibrum distribution function


    

In [4]:

    

#================ Initial boundary condition
w[0]=0.0
w[1]=0.5
w[2]=0.5  # k=3


    

In [5]:

    

#================== Initial value 
T[:]=0.0 #initial temperature
T[0]=100.0
for i in range (k): #k=0,1,2
 f[i,:]=w[i]*T[:]


    

In [6]:

    

#   Main loop  : comprised two parts :collision and streaming
#=====================
for n in range(mstep) :
# collision process
 for i in range(nx+1):
  T[i]=f[0,i]+f[1,i]+f[2,i]
  feq[0:k,i]=w[0:k]*T[i]
  f[0:k,i]=(1.-omega)*f[0:k,i]+omega*feq[0:k,i]
# streaming process
 for i in range(0,nx):
  f[1,nx-i]=f[1,nx-i-1]
  f[2,i]=f[2,i+1]

# Boundary condition 
 f[1,0]=Tleft-f[2,0]   # constant temperature at left  T=100
 f[:,nx]=f[:,nx-1] # adiabatic
 T[0]=100.0
 T[nx]=T[nx-1]
# end of the main loop


    

In [7]:

    

#     Finite difference # 
Tf=numpy.zeros(nx+1)   # Temperature finite difference
To=numpy.zeros(nx+1)   # Temperature storaage for previous time
Tf[0]=100.0

nx=30
dx=1./nx
dt=0.163  
nt=500              #  2*alpha*dt/dx^2  << 1 === > for stability
for n in range(nt) :
 Tf[nx]=Tf[nx-1]
 To[:]=Tf[:]
 for i in range(1,nx):
  Tf[i]=To[i]+ ( (dt*alpha/(dx**2.)) *(To[i+1]-2*To[i]+To[i-1]) )


    

In [9]:

    

plt.figure(figsize=(10,5), dpi=100)
plt.xlabel('x', fontsize=10) #x label
plt.ylabel('T', fontsize=10) #y label
plt.plot(x,T, 'r-',label=' Lattice Boltzmann Method'); 
plt.plot(x,Tf, 'bo',label=' Finite Difference Method '); 
plt.legend();


    

In [ ]: